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ABSTRACT 

We present a numerical study of the relation between the cosmic peculiar velocity field 
and the gravitational acceleration field. We show that on mildly non-linear scales (4-10 
h~ 1 Mpc Gaussian smoothing), the distribution of the Cartesian coordinates of each 
of these fields is well approximated by a Gaussian. In particular, their kurtoses and 
negentropies are small compared to those of the velocity divergence and density fields. 
We find that at these scales the relation between the velocity and gravity field follows 
linear theory to good accuracy. Specifically, the systematic errors in velocity-velocity 
comparisons due to assuming the linear model do not exceed 6% in (3. To correct for 
them, we test various nonlinear estimators of velocity from density. We show that a 
slight modification of the cv-formula proposed by Kudlicki et al. yields an estimator 
which is essentially unbiased and has a small variance. 
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1 INTRODUCTION 

According to the gravitational instability paradigm, struc- 
tures in the universe formed by the growth of small inho- 
mogeneities present in the early Universe. Gravitational in- 
stability gives rise to a coupling between the density and 
peculiar velocity fields on scales larger than the size of clus- 
ters of galaxies, the largest bound objects in the Universe. 
On very large, linear scales, the relation between the density 
contrast <5 and the peculiar velocity v in co-moving coordi- 
nates can be expressed in differential form, 



5(r) 



-(i/o/rV • v(r) . 



or in integral form, 
v(r) = g(r) . 
Here, 

dV 5{r') (r' - r) 



g(r) = H f 



-in 



(1) 



(2) 



(3) 



is a quantity proportional to the gravitational field, ex- 
pressed in units of km-s -1 . The coupling constant, /, carries 
information about the underlying cosmological model, and 
is related to the cosmological matter density parameter, Q m , 
and cosmological constant, Qa, by 



/(fim,fiA) 



70 



( 1 + t) 



(4) 



(Lahav et al. 1991). Hence, comparing the observed den- 
sity and velocity fields of galaxies allows one to constrain 
Q m , or the degenerate combination f3 = S7^ 6 /6 in the pres- 
ence of galaxy biasing (e.g., Strauss & Willick 1995 for a 
review). This comparison is done by extracting the density 
field from full-sky redshift surveys (such as the PSCz; Saun- 
ders et al. 2000), and comparing it to the observed velocity 
field from peculiar velocity surveys. The methods for doing 
this fall into two broad categories. One can use equation 
to calculate the predicted velocity field from a redshift sur- 
vey, and compare the result with the measured peculiar ve- 
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The linear relation between the peculiar velocity and the grav- 
ity, g, is (e.g., Peebles 1980) v = 2f/(3HSl m ) g. For simplicity, 
here we define the scaled gravity, g = 2//(3-fff2 m ) g. Then eq. J2J 
follows and g is given by eq. 13*1 , 
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locity field; this is referred to as a velocity-velocity compar- 
ison. Alternatively, one can use the differential form, equa- 
tion Q, and calculate the divergence of the observed veloc- 
ity field to compare directly with the density field from a 
redshift survey; this is called a density-density comparison. 
Nonlinear extensions of equation Q have been developed 
by a number of workers (Nusser et al. 1991, Bernardeau 
1992, Gramann 1993, Mancinelli et al. 1994, Mancinelli & 
Yahil 1995, Chodorowski 1997, Chodorowski & Lokas 1997, 
Chodorowski et al. 1998, Bernardeau et al. 1999, Dekel et 
al. 1999, Kudlicki et al. 2000, hereafter KCPR; see also the 
discussion below). However, very little work has been done 
to test on what scales the integral relation, equation 
holds, and how it might be extended into the mildly non- 
linear regime; thus the motivation for this paper. Attempts 
have been made to carry out velocity-velocity comparisons 
with very large smoothing lengths (e.g., the inverse Tully- 
Fisher method of Davis, Nusser, & Willick 1996), and very 
small smoothing lengths (the VELMOD method of Willick 
et al. 1997; Willick & Strauss 1998). Davis et al. (1991), and 
more recently Berlind et al. (2000) discuss the systematic 
errors caused by mismatch of smoothing scales between the 
velocity and density fields. In this paper, we concentrate on 
the velocity-velocity comparison after smoothing on scales 
of 4 /i _1 Mpc or larger: any smaller would be affected by 
strongly nonlinear effects, while larger smoothing would re- 
duce the number of independent volumes over which the 
comparison could be made. 

The amplitude of the velocity field smoothed on a given 
scale R with the window W depends on the density field 
power spectrum, P(k), as 



(v 2 ) oc 



dkP{k)W 2 {kR), 



(5) 



while for the density field the relation is as follows: 

{S 2 } oc / dkk 2 P(k)W 2 (kR) (6) 



(see the discussion in chapter 2 of Strauss & Willick 1995). 
Here, W is the Fourier transform of the smoothing window. 
The absence of the k 2 term in equation @ means that the 
velocity field is more heavily weighted by modes with low 
values of the wavenumber k, i.e., large scales which are fully 
in the linear regime. Therefore, we expect the relation be- 
tween the velocity, v, and the gravity, g, to be closer to linear 
than that between the density and velocity divergence. 
KCPR have shown that the relation between 



= -(H /) _1 V-v(r) 



(7) 



(see equation^ and 5 (proportional to V • g(r)) is nonlinear 
on small scales, and have proposed a semi-empirical formula 
accurately describing the dependence of 8g = (8\S) on 5: 

6s = a[(l + 5) 1/Q -l] +e. (8) 

Here, the constant e is introduced in order to keep (8) = 
as it must; it is approximately 



2a 
where cr| 



1 -*1 + O{of) 



(9) 



(5 2 ) denotes the variance of the density field. 



we use the exact value of e, obtained numerically. KCPR 
have found that a = 1.9 is a good fit over a large range of 
smoothing scales. Solving EquationQfor v in the case of an 
irrotational flow gives: 



v(r) 



Hof 
4n 



d 3 r , <Mr')(r'-r) 



(10) 



Thus, this approximation to e gets progressively worse as 
the fluctuations become more nonlinear and in the following 



where equation JHJ gives an expression for 8s. 

In this paper, we investigate the nonlinearities of the 
relationship between the velocity and gravity fields using a 
set of grid-based simulations. The simulations are described 
in §2. In §3, we ask how well the probability distribution 
functions (PDFs) of the Cartesian components of v and g 
are fitted with a Gaussian. Given that the source fields for 
the velocity and gravity fields (i.e. the velocity divergence 
and density contrast respectively) have mildly non-Gaussian 
distributions, it is not a priori obvious what the distribution 
of the integral quantities should be. In 0we directly mea- 
sure the relation between v and g on various scales, and test 
the extent to which linear theory, or non-linear extensions 
to it, may hold. This is important in determining whether 
the existing velocity-velocity comparisons which use linear 
theory give biased results. We present our conclusions in §5. 
Two appendices contain derivations of results used in the 
text. 



2 THE SIMULATIONS 

We performed our simulations using the CPPA (Cosmolog- 
ical Pressureless Parabolic Advection) code (see Kudlicki et 
al. 1996, KCPR). Matter in this code is represented as a 
non-relativistic pressureless fluid, and its equations of mo- 
tion are solved on a uniform grid fixed in comoving coordi- 
nates. Periodic boundary conditions are applied. Parabolic 
interpolation of the hydrodynamical state (like in the Piece- 
wise Parabolic Method scheme, see Colella & Woodward 
1984) assures low internal diffusion of the code and accurate 
treatment of high density and velocity gradients. 

We chose to use a grid-based code rather than an N- 
body code, because it produces a volume-weighted velocity 
field directly. This is important because equation 1101 is a 
solution to equation J7J only when v is a potential (curl- 
free) field, and the mass-weighted velocity field exhibits curl 
even in the linear regime. 2 Moreover, the a priori unknown 
galaxy bias does not allow one to treat observational data 
as purely mass-weighted anyway. 

All our simulations assume an Einstein-de Sitter uni- 
verse. The relation between the velocity and the (scaled) 
gravity in the mildly nonlinear regime is insensitive to the 
Cosmological Density Parameter and Cosmological Con- 
stant, as demonstrated both analytically (Gramann 1993, 
Chodorowski 1997, Nusser & Colberg 1998; see also App. B.3 



2 The mass-weighted velocity field is proportional to (1 + 5)v. If 
its curl is to be zero, with V X v = 0, one would require V<5 X v = 0. 
There is no a priori reason for this even in the linear regime. The 
mean value of the cosine of the angle between V<5 and v measured 
on the 30 h~ 1 Mpc scale in our six high-resolution simulations 
(see Table 1) is 0.81. The directions of these vectors are thus 
correlated, but not parallel. 
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Table 1. Parameters of the studied models. Hereafter, the model 
with the cell size equal to 1.56 h~ Mpc will be called the 'high- 
resolution model', while the remaining models will be called 'low- 
resolution models'. 



grid 


box size [h Mpc] 


cell size [h 1 Mpc] 


64 3 


200 


3.13 


128 3 


200 


1.56 


128 3 


400 


3.13 



of Scoccimarro et al. 1998), and by means of N-body simu- 
lations (Bernardeau et al. 1999), thus our results should be 
valid for any cosmology. 

Our simulations start from Gaussian density fluctua- 
tions with the linear APM power spectrum (Baugh & Efs- 
tathiou 1993, 1994, Baugh & Gaztafiaga 1996). The initial 
velocity field is obtained from equation The initial fields 
are evolved until the linear variance of the density in spheres 
of radius 8 /i~ 1 Mpc, <jg, is unity. Then for subsequent anal- 
ysis the output is selected for which nonlinear as = 0.87 
(the cluster normalization in the currently preferred model 
fl m = 0.3, Q.A = 0.7; Eke, Cole & Frenk 1996). 

In order to test the dependence of the results on the spa- 
tial resolution and box size we have studied three numerical 
models with parameters given in Table To improve the 
statistics, we have performed six realizations of each of the 
models, with different random phases of the initial density 
field. We have also performed a few additional simulations 
with the timestep reduced five times; the results did not 
change noticeably. 



1 1 I 1 1 

200 tr 1 Mpc, 64 3 
200 h" 1 Mpc, 12B 3 
400 tr 1 Mpc, 128 3 




Figure 1. The probability distribution of a Cartesian compo- 
nent of the 4 h~ 1 Mpc Gaussian-smoothed velocity field from our 
simulations and its deviations from Gaussianity. In both panels, 
the dotted line represents the results for the 64 3 model in the 200 
/i _1 Mpc box, dashed - the 128 3 model in the 200 ft. _1 Mpc box, 
and long-dashed - the 128 3 model in the 400 /i _1 Mpc box. The 
solid line shows a standardized Gaussian, G(v/a v ). 



3 THE MARGINAL DISTRIBUTIONS OF V 
AND G 

The typical correlation length of the density field is of the 
order of 5 /i _1 Mpc, while g is influenced by density fluctua- 
tions in a much larger region. For instance, the gravitational 
acceleration of the Local Group receives considerable con- 
tributions from shells up to at least 150 /i _1 Mpc (see e.g. 
Rowan- Robinson et al. 2000). Thus, g, and similarly v, come 
from integration over an effective domain containing a large 
number of essentially independent regions. Hence, the Cen- 
tral Limit Theorem suggests that they should have close to 
Gaussian distributions. 

We test the Gaussianity hypothesis with our simula- 
tions. We plot the measured distribution functions for in- 
dividual Cartesian components of the peculiar velocity and 
gravity fields, which we label v and g. It turns out that 
on mildly nonlinear scales, the distributions are well-fitted 
by Gaussians. Figures and |2] show the distributions for 4 
h~ Mpc Gaussian smoothing. (As there is no preferred di- 
rection in space, the distributions must be even functions. 
Therefore, without the loss of information we plot them as 
functions of the absolute value of a Cartesian component.) 
The closeness to a Gaussian distribution is remarkable, es- 
pecially for the velocity field. 

A standard way of quantifying modest departures from 
Gaussianity is to decompose the PDF with an Edgeworth 
expansion: a leading Gaussian, plus correction terms with 



amplitudes proportional to the higher-order connected mo- 
ments of the field (Longuet-Higgins 1963, Juszkiewicz et 
al. 1995, Bernardeau & Kofman 1995). As there is no pre- 
ferred direction in space, the skewness of the PDF of the 
Cartesian coordinates of the velocity and gravity fields must 
equal zero. Thus the first non-vanishing connected moment 
of a Cartesian component of the velocity field is its kurto- 
sis, K. v = {(v 4 ) — 3{v 2 ) 2 ) / (v 2 ) 2 , and similarly for the gravity 
field. Therefore, the kurtosis measures the leading-order de- 
parture from Gaussianity of the fields (Kofman et al. 1994, 
Catelan & Moscardini 1994). Specifically, the Edgeworth ex- 
pansion in our case reads: 



i + 



24 



(11) 



Here, P(fJ,) is the PDF for the variable v/a v or g/u g , where 
a v = {v 2 ) an d o 2 g — {g 2 } are the variances of the velocity and 
gravity fields respectively, and 4>(n) = (27r) -1 ^ 2 exp(— /x 2 /2) 
is the standardized normal distribution. The symbol Ha de- 
notes the fourth order Hermite polynomial. The quantity 54 
is related to the velocity, or gravity, kurtosis and the vari- 
ance of the density field in the following way: 



C 1 

04i> &s 



and 



■ Sigcr 2 . 



(12) 



The quantities S& v and iS4 S are called hierarchical ampli- 
tudes. For a given smoothing scale, during the weakly non- 
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Figure 2. The probability distribution of a Cartesian component 
of the 4 h~ 1 WLpc Gaussian-smoothed gravity field from our sim- 
ulations and its deviations from Gaussianity. Line coding is as in 
Figure Thick solid line in the lower panel shows the prediction 
of equation 1111 . with K, g = <S4 9 o"| = 0.83. 



Figure 4. The probability distribution of a Cartesian compo- 
nent of the 4 h~ -'-Mpc Gaussian-smoothed gravity field from high- 
resolution simulations, for it? = 0.25. In the upper panel, the 
dashed line shows the simulated PDF while the thin solid line - 
a standardized Gaussian. In the lower panel, the thick solid line 
shows the prediction of formula 1111 . with K, g = 0.22. 



linear evolution they are constant, independent of the nor- 
malization of the power spectrum. 

We plot K, v and K, g as a function of the Gaussian 
smoothing radius, R, in Figure[3] We find that IC V is close to 
zero for smoothing between 3 and 15 h~ Mpc. In contrast, 
the gravity field develops a detectable kurtosis. Still, this is 
substantially less kurtosis than the kurtosis of the gravity's 
source field, density. Measured from the high-resolution sim- 
ulations (i.e., with the grid size 128 3 and the box size of 200 
/i -1 Mpc), for the 4 ft - 1 Mpc smoothing the density kurtosis 
is ICs — 22.2 ± 5.0, thus more than an order of magnitude 
bigger than the gravity kurtosis. (For the same smoothing, 
the kurtosis of the velocity divergence is tCg = 4.1 ± 0.6). 
Therefore, qualitatively the picture is clear: our simulations 
reveal substantially less kurtosis of the velocity and gravity 
fields than those of the velocity divergence and density con- 
trast respectively; consequently the distributions of v and g 
are much closer to Gaussian. Our findings are inconsistent 
with the results of perturbative calculations of Catelan & 
Moscardini (1994) , and consistent with the results of N-body 
simulations of the velocity field of Kofman et al. (1994). 3 

3 The kurtosis of the velocity field, computed perturbatively by 
Catelan & Moscardini (1994), is greater than unity already for a 
Gaussian smoothing scale as large as 10 h~ Mpc for all cosmo- 
logical models they have considered, and grows with decreasing 
smoothing scale approximately like cr|(i?). In contrast, Kofman 
et al. (1994) find that 'the PDF of velocity (...) is almost indis- 
tinguishable from Gaussian in the simulations'. 



Quantitatively, the agreement of the evolved velocity 
and gravity PDFs with the perturbatively motivated equa- 
tion 111 It is less than perfect. In the lower panel of Figure |2] 
we plot the difference between the measured distribution 
for a Cartesian component of the gravity field and the stan- 
dardized Gaussian. As the thick solid line we plot the differ- 
ence predicted from equation 1111 . with the value K, g — 0.83 
(the mean value for all models) . The predicted difference has 
somewhat too high amplitude. 

Expansion 1111 is an expansion in the density variance. 
In our simulations, for 4 /i~ 1 Mpc smoothing erf = 0.65. 
For such a big variance, higher-order contributions to for- 
mula lllll may simply not be negligible. We have checked 
this conjecture by analyzing the gravity field at an earlier 
output time. In Figure we plot the gravity PDF for 
oj = 0.25. Formula lllll fits the simulated distribution bet- 
ter. 

We have also checked the perturbative scaling of the 
gravity (and velocity) kurtosis with the variance of the den- 
sity field, i.e. equation 1121 . Although for a given smoothing 
scale, during weakly non-linear evolution the hierarchical 
amplitude Si g remains constant, it depends moderately on 
the smoothing scale via an effective spectral index of the 
power spectrum at this scale. However, we have performed 
an additional simulation for a power-law initial spectrum 
(n = —1), and found qualitatively similar results; a change 
of sign in the kurtosis at several megaparsecs. 

Intrigued by this feature, we have studied the values of 
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Figure 3. Kurtosis of a Cartesian component of the velocity field (left) and gravity field (right) from our simulations as a function of 
the Gaussian smoothing length. The points are averages over six simulations times three directions, while the error bars are the standard 
deviations over this ensemble. 



the kurtoses for even larger, apparently linear, scales. The 
results for the gravity kurtosis are shown in Figure (The 
results for the velocity kurtosis are similar.) The values of 
the kurtosis clearly depend on the box size used in the sim- 
ulations. For the box size equal to 200 h~ 1 Mpc, they tend, 
instead to zero, to an asymptotic value —1.5. We interpret 
this as a purely numerical effect. Namely, on large scales the 
velocity field of the simulation is dominated by a small num- 
ber of Fourier waves, and thus the Central Limit Theorem 
does not force the distribution to be Gaussian, even in the 
initial conditions. Equations (HJ and © show that this ef- 
fect should be much stronger for g than for 8. On the largest 
scales, approaching the simulation box size, K, g converges to 
— 1.5, i.e. the value expected for a single mode, as shown 
in Appendix [X] 4 On scales shown in Figure simulations 
with the 400 ft -1 Mpc box size are much less affected by this 
effect, as expected. 

Why is it that this seems to affect scales down to ~ 20 
megaparsecs? Aren't there enough modes with scales rang- 
ing from ~ 20 megaparsecs to the box size, to warrant Gaus- 
sianity? The answer to this question is provided by equa- 
tion To the variance of the velocity field on a given 
scale contribute all larger modes with amplitudes propor- 
tional to P(k), the density power spectrum. The spectrum 
of the APM galaxies, employed by us, has a maximum at 



4 As pointed out by Scherrer (private communication), it is 
straightforward to show that the kurtosis of a Cartesian com- 
ponent of an isotropic gravity field satisfies Kg > —6/5, and sim- 
ilarly for the velocity kurtosis. Therefore, any observations which 
violate this bound indicate that one is in a regime in which the 
assumption of isotropy breaks down. This happens in our simula- 
tions with the box size of 200 h~ 1 Mpc for the smoothing scale of 
60 h~ 1 Mpc, at which K g < —1.2 (see Fig.|SJ. This is consistent 
with our idea that for such a big smoothing scale, a single mode 
dominates. 



200 tr« Mpc, 64 3 

200 h"» Mpc, 128 3 

400 h-< Mpc, 128 3 




20 



40 



60 



R (h-> Mpc) 



Figure 5. Kurtosis of a Cartesian component of the gravity field 
from our simulations as a function of the Gaussian smoothing 
length, now over a larger range of smoothing scales. The points 
are averages over six simulations times three directions, while the 
error bars are the standard deviations over this ensemble. Line 
coding is as in Figure^ At large scales, the results for simulations 
with the box size of 200 h~ 4 Mpc converge to the same values, 
regardless the resolution. 

~ 300 7l _1 Mpc. This is larger than the 200 fe _1 Mpc box 
size of some of our simulations; thus there are relatively few 
modes on the scale of the box size that dominate the velocity 
field on large scales. 

To settle this issue definitely would require additional 
simulations with box size comparable to the Hubble radius 
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(like the Hubble Volume simulations of the VIRGO consor- 
tium, Evrard et al. 2002) . 5 However, the ultimate goal of this 
paper is to study the relation between velocity and gravity 
for 4 h~ 1 Mpc smoothing. Figure [3] shows that on such a 
small scale the velocity and gravity kurtoses are positive, 
thus most likely induced by nonlinear dynamics rather than 
due to finite volume of our simulations. However, given the 
limitations of the Edgeworth expansion, 6 we would like to 
have an alternative measure of non-Gaussianity. As such a 
measure we borrow the concept of negentropy from informa- 
tion theory (e.g. Cover & Thomas 1991, Papoulis 1991). 
We define the entropy S[f] of a probability distribution 

S[f] = - f f(x)logf(x)dx. (13) 

It can be shown that, for a given variance, the entropy is 
maximal for Gaussian fields. Hence the difference between 
the entropy of a given field and the entropy of a Gaussian 
of the same variance, the negentropy Af[/], can be used as a 
measure of departure from Gaussianity: 

Af[/] = S[Gauss] - S[f] . (14) 

The negentropy was calculated by numerically integrating 
the integral of equation 1131 . using a PDF binned over the 
range from — 5a to 5a. We have found that this technique ap- 
plied to a Gaussian with the same range and binning yielded 
a negentropy of less than 10 -5 , which assures us that our 
results are not affected by numerical effects. 

We plot the results as a function of Gaussian smooth- 
ing scale in Figure The negentropy of the density field 
is compared to the negentropy of the gravity in the upper 
panel. The values for velocity and its divergence are com- 
pared in the lower panel. The negentropy of the velocity 
field is practically zero. The negentropy of the gravity field, 
even on the smallest scales, is extremely small compared to 
the negentropy of the density. The contrast between the ne- 
gentropy of the gravity and the negentropy of the density 
is even greater than between the corresponding kurtoses. It 
is so because the density field, unlike the gravity field, has 
significant skewness, also contributing to the value of the 
negentropy. In other words, the density field is more non- 
Gaussian than the gravity field not only because it has big- 
ger kurtosis, but also because it has non- vanishing skewness. 

We conclude that on mildly non-linear scales, the non- 
Gaussianity of v and g is completely negligible compared to 
the non-Gaussianity of 9 and 8. As stated earlier, this finding 
is consistent with the results of N-body simulations of the 
velocity field of Kofman et al. (1994). It is also consistent 
with the results of Gooding et al. (1992) and Scherrer (1992), 
where non-Gaussianity of the density field was (at least on 
large scales) due to initial conditions rather than nonlinear 
gravitational evolution. 



5 We have performed an additional simulation with the grid size 
256 3 and the box size of 800 h~ 1 yipc. On scales shown in Fig- 
ure 151 the gravity kurtosis tended asymptotically to zero, remain- 
ing positive. 

6 The Edgeworth expansion is an asymptotic expansion, not 
guaranteed to converge. It is not even guaranteed to be positive 
definite. 
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Figure 6. Negentropies of the cosmic fields, for the high- 
resolution simulations (i.e., with the grid 128 3 and 200 h~ 1 M.pc 
box size) . Upper panel compares the negentropies of gravity g and 
of density contrast 8, lower - negentropies of velocity v and of its 
divergence 9. 
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Figure 7. The joint probability distribution of a Cartesian com- 
ponent of the cosmic velocity and gravity fields, both smoothed 
with a 4 h~ 1 Mpc Gaussian filter. Solid contours represent com- 
bined numerical results from six high-resolution simulations, dot- 
ted ones - a standardized bivariate Gaussian of the same correla- 
tion coefficient. The contours correspond to the probability levels 
of 68%, 95%, and 99%. 



4 THE ONE-COMPONENT V-G RELATION 

In grid simulations, the shortest Fourier modes correspond 
to the Nyquist wavelength, of two cells. Smaller structures 
are not well resolved by the code. In low-resolution simu- 
lations (cell size 3.13 ft _1 Mpc, see Table 0, the Nyquist 
wavelength is larger than the scale at which we would 
like to study the velocity-gravity relation (4 /i _1 Mpc). In 
high-resolution simulations (cell size 1.56 /i _1 Mpc), it is 
smaller. Therefore, to model the velocity-gravity relation 
at 4 /z -1 Mpc scale, we will use only high-resolution simula- 
tions. 

On scales of a few /i _1 Mpc, the relation between 8 and 
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6 is non-linear. Nevertheless, since the probability distri- 
butions of the Cartesian components of both v and g are 
nearly Gaussian, we hypothesize that their joint probability 
distribution is a (bivariate) Gaussian as well. In Figurc[7|wc 
present the simulated joint PDF for v and g, measured on 
the 4 h~ Mpc scale. It is indeed quite close to a bivariate 
Gaussian. The only substantial deviation is for the proba- 
bility contour of 99%, along the gravity coordinate. There, 
positive kurtosis of the gravity field broadens the simulated 
isocontour with respect to the Gaussian one. Figure |7| also 
shows that v and g are strongly correlated. In Add. IB1 we 
show that in the bivariate Gaussian case, the relationship 
between v and g is purely linear; we now show that this is 
approximately the case in the simulations. 

In the linear regime v = g and thus each of the Carte- 
sian components of these quantities, which we denote v and 
g, respectively, are also equal. In practice, the relationship 
between these two quantities has some finite scatter (Fig.0, 
thus we will characterize the relation between the mean v at 
constant g, (v\g), and the converse, (g\v). Fitting a straight 
line to (v\g) as a function of g gives a slope of 0.94. There- 
fore assuming pure linear theory on this scale would give a 
systematic 6% error in /3. We now consider going beyond 
linear theory. 

The relationship between (v\g) and g should be invari- 
ant to coordinate inversions, thus it must be an odd function; 
similarly for the relationship between (g\v) and v. 

Hence, we shall adopt third-order polynomials 



1.01 



(g\v) = civ + C3V 3 /al 
(v\g) = dig + d 3 g 3 /cTg 



(15) 
(16) 



as the simplest odd non-linear model, and use the unitless 
parameters C3 and dz as a measure of the non-linearity of 
the relation. As the deviations from linear theory are small, 
we expect ci and d\ to be close to unity, and C3 and d$ 
to be small. Since in velocity-velocity comparisons one re- 
constructs velocities based on the gravity field, the quan- 
tity (v\g) is more relevant, and we shall concentrate on it 
in this paper. Note however that relation 115H can be used 
to transform velocities in the first step of density-density 
comparisons such as POTENT (Dekel et al. 1999). 

The parameters ci and C3 as functions of the smoothing 
scale are shown in Figure [S] while Figure [5] shows di and d$. 
We have fitted them independently for the three coordinates 
in each of the six realizations of every numerical model; the 
lines are the averages over these 18 datasets and errorbars 
indicate their standard deviations. As expected, C3 and d$ 
are much smaller than unity. On large scales, ci and di tend 
to the linear theory value with accuracy better than half per 
cent. Moreover, they remain close to unity even on small 
scales, which implies that the systematic bias in estimating 
(3 based on the linear theory approximation is small. 

To quantify this bias, in Figure HUI we show the param- 
eters ci and di when C3 and ^3 are set to zero. For R > 20 
/i _1 Mpc, c\ and di deviate from unity by less than 2%. 
Therefore, at these scales one can apply linear theory with 
such good accuracy. 

Note in Figures |H| and [5] that at small scales, high- 
resolution simulations yield the highest values of C3 and 1— ci 
(similarly for |e?3| and 1 — di), but overall the effects of res- 
olution and of the simulation box size are not large. 
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Figure 8. The parameters c\ and C3 of the polynomial approxi- 
mation to the mean one-component v— g relation (Eq. 1151 as func- 
tions of the smoothing scale. Line coding is as in Figure HI 
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Figure 9. The parameters di and ds of the polynomial ap- 
proximation to the mean one-component g— v relation (Eq. 1161 as 
functions of the smoothing scale. Line coding is as in Figure HI 



We plot the mean predicted velocity as a function of the 
true velocity for a 4 /i _1 Mpc Gaussian filter in Figure FiTI 
The velocity predicted with the linear-theory model (i.e., the 
gravity) is plotted as a dashed line. Combined data from six 
realizations of our high-resolution model (three components 
each) are binned with respect to the predicted velocity and 
averaged in the bins. All velocities are in kilometers per 
second (the one-dimensional rms true velocity a v = 290 km- 
s _1 on this scale). A long-dashed line shows the velocity 
predicted by the polynomial formula 1161 . We see that linear 
theory, v = g, is a good fit for velocities up to about 2a v . The 
polynomial approximation works well up to about 3cr„, but 
fails in the tails too. Therefore, we decided to seek a better 
estimator of velocity from density than a simple function of 
gravity. 
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Figure 11. One Cartesian component of the velocities recon- 
structed from the density field - combined data from six high- 
resolution simulations, observed with a 4 h~ 1 Mpc Gaussian filter. 
Top: mean predicted velocity given true velocity. Bottom: errors 
of the reconstruction. All velocities are in kilometers per second. 
The one-dimensional rms true velocity, cr v , for this smoothing 
equals 290 km-s -1 . Dotted line: identity, dashed line: linear the- 
ory model, long-dashed: a polynomial model (Eg. 1161 . Solid line 
shows the reconstruction from the density field computed using 
the non-linear 07- formula (Eq.E3, with a = 1.56 and 7 = 1.054. 



KCPR have found that formula @ is a good description 
of the tails of the relation between the density and the ve- 
locity divergence, so it should also work well on the integral 
level feci. HUP . It turned out, however, that at small scales 
a small modification is needed. Formula © implies that 
the coefficient of the linear term in density in the density- 
velocity divergence relation is unity, 6s = S— r2(S 2 — erf ) + . . ., 
where r-z = (a — 1)/2q, regardless the value of a, while we 
have found here that at small scales it is not strictly true. To 
cure this problem, we introduced an additional parameter, 
7, so that 

-5) 1/c 



= 07 [(1 



-1 +e. 



(17) 



For given a and 7 the predicted velocity, u pr od, was cal- 
culated from equation 1101 . The best values of a and 7 
were found by minimizing the quantity X 2 = J^iC^pred — 
v t l ue ) 2 /(3N 3 ), where i>tr UC is the true velocity and the sum 
is over all grid points (iV 3 in total) and all Cartesian com- 
ponents. The resulting values of the parameters, for the 
smoothing scale 4 /i _1 Mpc and different values of as at dif- 
ferent output times, are shown in Figure IT21 

In the limit of linear theory, = 7 = 1. The offset 
e, fully determined by a and 7, is then equal to zero. In 
general, on the integral level the value of e is irrelevant, since 
its contribution to velocity in the integral in equation 1 101 
averages out to zero. We see that for small as the parameter 
7 tends to unity, while a decreases only weakly and remains 
well away from this value. This is consistent with the findings 
of KCPR, that a-formula (eq. |H] or eq. El with 7 = 1) 
describes well the density-velocity relation in the weakly 
nonlinear regime (i.e., for erg smaller than, say, 0.3). Though 
the best-fit value of 7 is close to unity even for large as, its 
inclusion markedly decreases X 2 . 

In Figure EH we plot the quantity x(v) = \/AX 2 /AiV, 
that is the square root of the contribution to X 2 from a given 
bin in velocity, divided by the number of points in that bin. 
(X has units of km • s _1 .) We have X 2 = J \ 2 ( v ) P( v ) iv, 
where P(v) is the probability distribution function of a 
Cartesian component of the velocity field. Thus, X 2 is a 
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Figure 13. The quantity \ (see text). Dashed line shows it for 
the velocity predicted according to the linear theory, long-dashed 
— for the velocity approximated by a polynomial in the gravity, 
and solid - for the velocity predicted by the «7-formula (inserting 
equation 1171 in I1UI . 



number- weighted average of \ 2 ■ Dashed line shows % f° r 
the velocity predicted according to the linear theory, long- 
dashed - for the velocity approximated by a polynomial in 
the gravity, and solid - for the velocity predicted by the 07- 
formula (inserting eg. H7l in llUII . The linear-theory estimator 
of velocity yields the largest values of \- The polynomial 
formula results in smaller \ m the tails, but in equal \ f° r 
'typical' values of velocity (— a v < v < a v ). The Q7-formula 
yields the smallest x m the whole range of velocity. For typ- 
ical values of velocity, i.e. where most of the data comes 
from, x is slightly smaller than for the previous estimators. 
Furthermore, it grows only moderately in the tails. 

In Figure fTTI a solid line shows the mean velocity pre- 
dicted by the Q7-formula as a function of the true velocity. 
The smoothing scale is that of velmod, i.e. 4 /i _1 Mpc, and 
erg = 0.87. For this value of as, the best-fitted values of a 
and 7 are a = 1.56 and 7 = 1.054. We see that the 07- 
formula yields a practically unbiased estimator of velocity. 
The difference between the predicted and true velocity in 
the tails is very small. 

Figure 1141 demonstrates that the estimator of velocity 
from density provided by the ct7-formula is not only essen- 
tially unbiased but also has a small variance (as suggested 
already by Figure 1131 . The figure shows a joint PDF for 
the velocity predicted by the Q7-formula and the true ve- 
locity, for 4 h~ 1 Mpc smoothing. It is instructive to compare 
it with Figure [7] where i> pre d = g- We have found that for 
larger values of the smoothing scale the correlation between 
the predicted and true velocity is even higher. 

In Figure 1151 we show the dependence of a and 7 on 
the smoothing scale, for as = 0.87. For large R the param- 
eter 7 tends to unity, as expected. The parameter a grows 
from the value around 1.55 for 3 /i _1 Mpc, to around 2.05 for 
20 /i _1 Mpc. This is in contrast with the results of KCPR, 
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Figure 14. The joint probability distribution of a Cartesian 
component of the velocity predicted by the o;7-formula and the 
true velocity, both smoothed with a 4 ft _1 Mpc Gaussian fil- 
ter. Solid contours represent combined numerical results from six 
high-resolution simulations, dotted ones - a bivariate Gaussian 
of the same correlation coefficient and variances. The contours 
correspond to the probability levels of 68%, 95%, and 99%. 
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Figure 15. The dependence of a and 7 on the smoothing scale, 
for as = 0.87. 

where a ~ 1.9 was a good fit for a large range of scales. We 
can see at least three likely sources of this discrepancy. First, 
the power spectrum used by KCPR was a pure power law, 
P(k) oc AT 1 , while ours is not and the effective spectral index 
of our APM spectrum depends (slightly) on the smoothing 
scale. Secondly, KCPR studied the a-formula, equation |H| 
which is equivalent to our 07- formula only when 7=1. Fi- 
nally, KCPR's fit was obtained for all bins in velocity hav- 
ing equal weight, while we weight the bins by the number of 
points. As a result, our fit is much less affected by the tails 
of the velocity distribution. 

To summarize, we have tested various estimators of ve- 
locity from density. Nonlinear effects can be accounted for ei- 
ther at the integral level, employing equation llul , or at the 
differential level, employing equation (11711 . We have shown 
that the differential method works much better. 
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5 SUMMARY AND CONCLUSIONS 

In this paper we have studied the cosmic velocity-gravity 
relation. First, we have measured the non-Gaussianities of 
the cosmic velocity and gravity fields, evolved from Gaussian 
initial conditions, by computing their kurtoses and negen- 
tropies. We have shown that, on scales of a few /i _1 Mpc, 
the non-Gaussianities of the cosmic velocity and gravity 
fields are small compared to the non-Gaussianities of ve- 
locity divergence and density. A similar finding for the ve- 
locity field was reported by Kofman et al. (1994). Guided by 
this result, we have shown that the relation between v and 
g is nearly linear. Moreover, its proportionality coefficient 
is close to that predicted by linear theory. Specifically, we 
have shown that the systematic errors in velocity-velocity 
comparisons due to assuming linear theory do not exceed 
6% in /3. (Strictly speaking, if erg < 0.9 and the smooth- 
ing scale is not smaller than 4 h~ 1 Mpc). To correct for this 
small bias, we have tested various nonlinear estimators of 
velocity from density. We have shown that the Q7-formula 
(a slight modification of the a-formula proposed by KCPR) 
yields an estimator which is essentially unbiased and of small 
variance. 

The smoothing of observed peculiar velocity data, with 
its sparse and noisy coverage of the velocity, is techni- 
cally difficult. Thus velocity-velocity analyses like velmod 
compare unsmoothed peculiar velocities with minimally 
smoothed predicted velocity fields from redshift surveys. 
The finite resolution of a grid code does not allow us to 
test this effect; it would be worthwhile to repeat the calcu- 
lations with a high resolution iV-body code, or with CPPA 
enhanced with Adaptive Mesh Refinement. The former has 
been done by Berlind et al. (2000), but they did not sep- 
arate the effects of non-linear evolution from the effects of 
different smoothing of the velocity and gravity fields. 

An alternative method of estimating from cosmic ve- 
locities is by comparing the cosmic microwave background 
dipole with the density dipole inferred from a galaxy redshift 
survey. In maximum-likelihood approaches to this problem 
(Strauss et al. 1992, Schmoldt et al. 1999, Chodorowski & 
Ciecielag 2002) the relevant quantity is the joint distribution 
for the velocity and gravity, which is commonly approxi- 
mated as a multivariate Gaussian. We have shown that this 
indeed holds to fairly good accuracy (Fig. |7J . Small nonlin- 
ear effects can be corrected for by a more thorough mod- 
elling of the joint distribution. There is, however, a better 
approach: from a redshift survey, instead of calculating the 
acceleration on the Local Group (i.e., the density dipole, or 
gravity), one can calculate the predicted velocity from the 
a7-formula. The joint distribution for the predicted veloc- 
ity and the true velocity is perfectly Gaussian (see Fig. lMH . 
Moreover, the relation between these two variables is tighter 
than the relation between the velocity and the gravity. This 
implies that in the proposed method, the estimated value of 
f3 will be unbiased and its inferred errors will be smaller. 
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APPENDIX A: CALCULATION OF THE 
KURTOSIS OF VELOCITY FOR A SINGLE 
FOURIER MODE 

Let us assume a single-mode distribution of the density field. 
All variables will depend on one dimension only. In the linear 
regime the density field is 



5{x, y, z) = S(x) — sin(a;) . 



(Al) 



(The amplitude of the wave is irrelevant here, since it will 
cancel out in the calculation of the dimensionless kurtosis.) 
Then, in the linear regime, 



8 = = -V- v = . 

dx 



(A2) 



[8 denotes the scaled velocity divergence, so in equation 1A2I 
there is no /(Sl m ,fiA) term.] Hence, 



,=-j S{x)dx 
and 

v(x) = cos(a:) . 



(A3) 



(A4) 



Let P(v) denote the volume-weighted probability distribu- 
tion function of a single component of the velocity field: 



P(v)\dv\ = -\dx\ . 

7T 

Hence, 
P(v) = 



7rsin(a;) _ cos 2 (x) 

and finally 

p(v) = TV- 1 (i - v 2 y 1/2 ■ 



(A5) 
(A6) 

(A7) 



Now let us compute the four first moments of this distribu- 
tion: 



(v) = a 
I 2 \ 1 

(v ) 



1 2 , 

v dv 



J-i Vl-v 2 



(v 3 ) = 

(V > = - 



v dv 



vr J_ 1 s/l - v 2 8' 
We can calculate the kurtosis of this field now: 
(v*) - 3{v 2 ) 2 3 



K v 



(A8) 
(A9) 
(A10) 
(All) 

(A12) 



= exp 



(p 2 — Irptv + v 2 ) 



2(1 -r 2 



(B2) 



Were the fields uncorrelated (r = 0), p{fi,v) would be just 
a product of two Gaussian distributions of fi and v. 

Since the variables are assumed to be correlated but in 
general not identical, the relation between them will have 
a scatter. The mean relation between [i and v is defined as 
mean fi given v, {fJ,\v). By definition, the conditional PDF 
of fi given u, p(fj,\v), is p(fi,v)/p(v). Equation 1B2I yields 



p(n\v) = 



yj2n(l - r 2 ) 



exp 



2(1 -r 2 



(B3) 



We see that the conditional PDF is a Gaussian with modified 
mean and variance. Specifically, 



(p\v) = ru, 
and 

{p?\v) - (fi\u) 2 = l-r 2 



(B4) 



(B5) 



The relation between two correlated Gaussian variables 
is thus linear. If the variables are uncorrelated there is no 
relation whatsoever. 7 Therefore, the relation between two 
Gaussian variables, if it exists at all, is always linear. 

One may define an 'inverse' relation to that specified in 
equation HB4II . i.e., (v\n). From symmetry of the joint PDF, 



(u\n) = r/i. 



(B6) 



The proportionality coefficient in the 'inverse' relation is 
thus not a simple reciprocal of the coefficient in the 'for- 
ward' relation, equation HB41 . The difference is equal to 



(1 — r )/r — 1 — r for r close to unity. In other words, 
the difference between the true and the 'naive' coefficient of 
the inverse relation is directly related to the scatter in the 
relation. 



APPENDIX B: THE MEAN RELATION 
BETWEEN TWO GAUSSIAN VARIABLES 

Consider two correlated Gaussian variables, yi and y2, with 
zero means and variances a\ = (y 2 ) and a 2 = (y 2 ), respec- 
tively. Their correlation coefficient is 



0"lO"2 



(Bl) 



Let us introduce normalized variables, [i = yi/ai and 
v = 2/2/02. We require that their joint probability distribu- 
tion function (PDF) is a bivariate Gaussian, 



7 This is only the case for Gaussian variables. For a counterex- 
ample in the opposite case, consider y\ and 1/2 = y 2 , where j/i is 
Gaussian. 
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